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Abstract: We have made a Markov Chain Monte Carlo (MCMC) analysis of primordial 
non-Gaussianity (/nl) using the WMAP bispectrum and power spectrum. In our analysis, 
we have simultaneously constrained /nl and cosmological parameters so that the uncer- 
tainties of cosmological parameters can properly propagate into the /nl estimation. Inves- 
tigating the parameter likelihoods deduced from MCMC samples, we find slight deviation 
from Gaussian shape, which makes a Fisher matrix estimation less accurate. Therefore, 
we have estimated the confidence interval of /nl by exploring the parameter likelihood 
without using the Fisher matrix. We find that the best-fit values of our analysis make a 
good agreement with other results, but the confidence interval is slightly different. 
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1. Introduction 

Over the past years, there have been great successes in measurement of the CMB anisotropy 
by ground and satellite observations [1-5]. CMB anisotropy, which is associated with 
the inhomogeneity of the last scattering surface, provides invaluable information on our 
Universe. In particular, the angular power spectrum of CMB anisotropy provides strong 
constraints on cosmological models and cosmological parameters [6-8]. 

Large classes of cosmological inflation models predict some level of primordial non- 
Gaussianity (/nl) ; whose detailed shape and amplitude depend on the particular class of 
the inflationary model [9-11]. Therefore, an investigation of the imprints of primordial non- 
Gaussianity allows us to constrain and rule out some classes of inflationary models [10-12]. 
Since the release of the COBE and subsequently WMAP data [1, 13-17], there have been 
active investigations on the imprints of primordial non-Gaussianity, using CMB data from 
the satellite and sub-orbital experiments [11,16,18-29]. Currently, the strongest limit is 
imposed on the local-type / NL = 32 ± 21 by the WMAP data [11,28]. CMB bispectrum, 
which is most sensitive to /nl parameters, also have some dependence on cosmological pa- 
rameters [9,11,19]. Noting this, we have simultaneously constrained /nl and cosmological 
parameters by using the WMAP bispectrum and power spectrum. In this way, we have 
allowed the uncertainties of cosmological parameters to properly propagate to the /nl es- 
timation. The parameter likelihoods deduced from MCMC samples show slight deviation 
from Gaussian shape, which makes Fisher matrix estimation less accurate. Therefore, we 
have estimated the confidence intervals of /nl by fully exploring the parameter likelihood 
without using Fisher matrix. We find our estimation makes a good agreement with the 
existing results, but the confidence interval of /nl is slightly different. 
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The outline of this paper is as follows. In Section ||, we discuss the primordial non- 
Gaussianity. In Section || and |], we discuss the CMB bispectrum and its likelihood function. 
In Section |5[ we make an analysis of WMAP data and present the results. In Section |6|, 
we summarize our analysis and discuss prospects. 



2. Primordial perturbations 

Primordial perturbation in Fourier space is given by [9,22,30,31]: 

$(k) = $ L (k) + $ NL (k), 



(2.1) 



where ^(k) and <l?NL(k) denote the primordial perturbation of purely Gaussian nature and 
the deviation respectively. In most of inflationary models, <&L(k), which follows a Gaussian 
distribution and statistical isotropy, have the following statistical properties [6-8]: 



<$ L (ki)$ L (k 2 )> = (27r) 3 P$(fc)<f(ki+k 2 ) 



(2.2) 



where (. . .) denotes the average over an ensemble of universes, and P$(k) = Ak n ~ is a 
primordial power spectrum with a normalization amplitude A. The non-Gaussian part of 
the primordial perturbation, $NL(k), has negligible 2-point correlation [12,32], which is 
also severely constrained by observational data: 



($NL(ki)$NL(k 2 )> « 0. 



(2.3) 



Additionally, we may consider a 3-point correlation. For Gaussian <l?L(k), we expect a 
vanishing 3-point correlation [6-9]: 



($ L (k 1 )$ L (k 2 )$ L (k 3 )) =0. 



(2.4) 



On the other hand, a non- vanishing 3-point correlation of <J>NL,(k) is predicted by large 
classes of inflationary models: 



($ NL (ki) $ NL (k 2 ) $N L (k 3 )> = (27r) 3 5(k! + k 2 + k 3 ) F{h, k 2 , k, 



3 ■ 



(2.5) 



where F(ki, A; 2 , /c 3 ) denotes a primordial bispectrum, whose detailed shape depends on the 
specific inflationary model [10]. In this paper, we shall consider the "local" and "equi- 
lateral" forms of the bispectrum, which are respectively associated with the squeezed and 
equilateral configuration of three vectors [9-11,22,28,30,31,33,34]. The primordial bispec- 
trum in its "local" and "equilateral" forms is defined as follows [10,11,35]: 
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where /^l anc ^ /ni!" denote the non-linear coupling parameters for "local" and "equi- 
lateral" forms respectively. The "local" form possesses a special significance because all 
single-field inflation models predict a negligible f^? al - Therefore, the detection of a high 
/nl*' wou ld allow us immediately to rule out single-field inflation models [11,36-39]. Be- 
sides the aforementioned forms, there are also other forms for primordial bispectrum, which 
are predicted by certain inflation models [40-43]. Therefore, the physical observables asso- 
ciated with the 3-point correlation provides important information on our early Universe. 



3. The CMB anisotropy 

The CMB anisotropy over a whole-sky is conveniently decomposed in terms of spherical 
harmonics: 



T(n) = y~]g lm Yi m (ti), 



Irn 

where a/ m and Y; m (k) are a decomposition coefficient and a spherical harmonic function. 
The decomposition coefficient a/ m is related to the primordial perturbations ^(k) according 
to: 

d 3 k 



$(k)«^)^(k), (3.1) 



where gi(k) is a radiation transfer function. Largely, we may consider CMB power spectrum 
and bispectrum, which are associated with 2-point and 3-point correlation of primordial 
perturbation respectively. In the absence of tensor perturbation, the expectation value of 
CMB power spectrum is given by: 

Q = {ai m a* m } (3.2) 
= ~y k 2 dkgf(k)P^(k). (3.3) 

As shown in Eq. |2.3| , 2-point correlation of $NL(k) is insignificant. Therefore, the CMB 
power spectrum is mainly associated with <3?L(k) [12,32] (c.f. Eq. |2.2j ). 

Considering the rotational properties of spherical harmonics, we can construct a rota- 
tionally invariant estimator of the CMB bispectrum as follows [44]: 

B ?ihh = y2 ( 1 2 3 ) a hm l ai2m 2 ai 3 m 3 . (3.4) 

^ t \ m m m j 

Hereafter, the condition (Zi < I2 < h) will be implied, unless stated otherwise. The 
expectation value of the CMB bispectrum, in other words, the theoretical CMB bispectrum, 
can be split into quantities associated with the different possible shapes of primordial non- 
Gaussianity: 



/ nobs \ \ fi r>i 
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where i denotes "local" and "equilateral" in our analysis. The CMB bispectrum of f^ L = 1 
is given by [9,11,22,28,34]: 



oo 

BjS = 2I hhh I r 2 dr [a h (r)(3 h (r)(3 h (r) + 2 perm.] , (3.5) 
Jo 

POD 

Khh = ~ W hhi + 6 4fafa / r 2 dr \p b (r) 7/2 (r)S h (r) + (5 perm.) - 26 h (r)S h (r)5 h (r)] , 

« 

(3.6) 

where "perm." denotes cyclic permutation and 



(3.8) 
(3.9) 
(3.10) 
(3.11) 
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with ji(x) being a spherical Bessel function. In our analysis, we have computed Eq. 3.5 and 



3.6 by numerical integration. We have investigated the functional shape of the integrand 
and used 10, 10, 100 and 30 integration points for the intervals r < 6.5, 6.5 < r < 13, 
13 < r < 16.5 and 16.5 < r < 39.5 [Gpc] respectively. Compared with a much denser 
sampling (x50), we find that the difference in the estimation of f}§f is less than 1%. 
For k space integration in Eq. 3^, |3.9| , 3.10 and p. 11 , we used the k space sampling of 



CAMB [45]. As discussed in Eq. 2.4 , 3-point correlation of ^(k) is expected to be zero. 
Therefore, CMB bispectrum is not associated with $L(k), but only with $NL(k) (c.f. Eq. 
P). 

Due to foreground contamination, CMB anisotropy cannot be measured reliably over 
the whole sky [17,46,47]. Given the spherical harmonic coefficients of a cut sky, we can 
construct the following estimator: 

Sobs _ ( h h h \ ~ ~ ~ 

Z A 1 m m m 1 

mim2m 3 \ / 

where 

dim =Jdn Y? m (n) T(n) W(n), (3.12) 

with W(n) being a foreground mask function. For the temperature-only case, the cut sky 
induces mostly a monopole contribution outside the mask, which can be made small by 
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subtracting the monopole term outside the mask [35,48]. To a good approximation, the 
CMB bispectrum estimated from a cut sky is related to that of a whole sky as follows: [9]: 



/ Bobs \ ^ f 1 nobs \ 

where / s k y is the sky fraction. 

In the weak non-Gaussian limit, the covariance matrix of Bf^i 3 is nearly diagonal and 
its diagonal elements are [11]: 

((AB$k) 2 > « fly( C h + N h )(C h + N l2 ){C h + N l3 )A hl2l3 (3.13) 

where C/ and Ni are the CMB and noise power spectrum, and Ay 2 / 3 is 1, 2 and 6 for all 
different three li, two equal l{ and l\ = I2 = h respectively. In the WMAP noise model, 
the pixel covariance is given by: 

N « = ^) s «' < 3 - 14) 

where gq is the rms noise per observation, and N Q ^ s (h) is the number of observations for 
the ith pixel of the sky direction n. Under a white noise assumption, one can estimate 
noise power spectrum Ni as follows [22]: 

n.^f«pm, ( 3. 15 , 

J 4vr/ sky iV obs (n) 



where Q pix is the solid angle per pixel. In deriving Eq. [3.13 , we have assumed that the 



spherical harmonic coefficient of the instrument noise rji m has a diagonal covariance (i.e. 
( r ]lm'ni> m /) = ^U'S mm 'Ni), which is true only for the case of homogeneous noise and a whole 
sky coverage. However, Monte-Carlo simulations show that non-diagonal parts of (fjimVi'm') 




Figure 1: Noise correlation (rjimVim') I {VlmVim) ^ or ' = ^00: real part (left) and imaginary 
part(right), estimated from 5 x 10 5 simulation of inhomogeneous noise on a cut sky. 

tend to be less than 1% of the diagonal parts. To be specific, we have simulated 5 x 10 5 
noise maps for the V and W band respectively by the WMAP noise model (c.f. Eq. 3.14| ), 
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and applied the foreground mask KQ75 to them. From the cut-sky simulations, we have 
estimated the noise covariance matrix {f]i m rjti f ) . The magnitude of noise correlation turns 
out to decrease with the relative distance in spherical harmonic numbers (i.e. \l — l'\ and 
\m — m'\). In Fig. [l], we show the histogram of noise correlation for the fixed multipole 
I = V = 400 and m ^ m' . Note that the values of the y axis denote the number of m and 
m' pairs, whose correlations correspond to the values of the x axis. As shown in Fig. [l], 
we find that non-diagonal parts of {i]im r li'm') tend to be less than 1% of the diagonal parts. 
Besides, we shall be estimating the bispectrum in the multipole range I < 400, where the 
Signal-to-Noise Ratio is high (e.g. 17 < Q/Ni < 42 at I < 400) [49]. Therefore, we find 



Eq. |3.13| provides good approximation for the variance. 

Point sources, whose fluxes are below a certain detection threshold, are not properly 
masked out, and therefore require an accurate assessment of their contribution to the 
bispectrum. Assuming a Poisson distribution for the point sources [50], their contribution 
is given by [9,11,22]: 

Kclhhh, (3-16) 
where 6 S rc is an amplitude to be determined, and Ii x i 2 i 3 is given by Eq. |3.7[ 

4. Likelihood function 

As shown in Eq. |3.4| , CMB bispectrum observables B°^i 3 contain a sum of a large number 
of terms. In the weak non-Gaussian limit, the correlation between the terms is insignificant: 



. / l\ 12 I3 \ I l\ 12 k \ I 

( / / / ] a hm 1 ai 2 m2 a l 3 m 3 (%m! a l 2 m' 9 a l 3 m'J } 

\ mi m 2 ^3 / \ IT^l 1^2 m 3 J 

(C h + N h )(C h + N h ){C h + N h )A hhl3 



0, (4.1) 



where we referred to Eq. |3.13j . Accordingly, the distribution of B°^f ls , which contains a 
large number of weakly correlated random variables, tends toward a Gaussian function by 
the central limit theorem. It is worth to note that, unlike CMB power spectrum observables, 
bispectrum observables B°^i 3 contain a sufficiently large number of degree of freedom 
(e.g. 125) even for the lowest multipole (e.g. 1 = 2). Therefore, we may approximate the 
likelihood function of the CMB bispectrum as follows: 

C{B hl2h {\ a )\B^ h ) = 1 r exp[-i(B obs - B(A a ))t C" 1 (B obs - B(A Q ))], (4.2) 

(27r) 2- |C| 2 z 

where C is a n x n covariance matrix, X a are the parameters: {/nl 3 " 1 ; /ni!" 1 ' ^src} + cos- 
mological parameters, and B obs and B(A a ) are vectors consisting of the observed and the 
theoretical bispectrum respectively. The theoretical bispectrum B(A Q ) is given by: 

B hhh {X a ) = w h w h w h /^[/^B£# 8 (A«) + fT l B^ 3 + b SIC I hhl3 ), (4.3) 

where wi is the effective window function, which is the product of the beam transfer function 
and the pixel window function. 
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Figure 2: Correlation between CMB power spectrum and bispectrum: (Bi lt i 2t i 3 Ci) 



In order to derive a total likelihood, we have investigated the correlation between 
CMB power spectrum and bispectrum. Due to the complexity of analytic derivation, we 
have made numerical investigation, using simulations produced by [51]. The simulations 
are provided in such a way that we may set the value of f§jf. In Fig. §, we show the 
correlation estimated from 1000 simulations, where /nl^ * s se ^ to 40. As shown in Fig. [2], 
we find the correlation between the CMB power spectrum and bispectrum is insignificant. 
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Therefore, to a good approximation, we may construct a full likelihood function as follows: 



£(Q\C? hs ) x C(B hhh (X a )\Bfl 



obs ^ 



(4.4) 



where C{Ci\C{) is the likelihood function associated with CMB power spectrum, whose 
approximate and optimal expression is found in [52]. The corresponding log-likelihood is 
given by: 



log C(Q | C? hs ) + log £{B hl2h (K)K? 2h ) 



(4.5) 



The log-likelihood of the WMAP CMB power spectrum, which corresponds to the first 



term in Eq. 4.5, is provided by the WMAP team's code and well-integrated into the widely 
used Monte-Carlo sampling code (i.e. CosmoMC). Therefore, we are going to be mainly 



concerned with adding the second term of Eq. 4.5 to the existing code. 



As shown in Eq. 4.5, the likelihood function has dependency on several cosmolog- 



ical parameters. Therefore, we are going to explore multi-dimensional parameter space 
via MCMC sampling. In this way, we are going to simultaneously fit all cosmological 
parameters, including /nl- 



5. Analysis of the WMAP data 

We have co-added WMAP 7 year foreground-reduced maps of the V band and the W 
band, and applied the KQ75 mask to the co-added map [1,46,47]. From the masked 
V+W map, we have estimated B°^i 3 up to / < 400 via Eq. |3.4| , where the Wigner 3j 
symbols have been computed efficiently using the recurrence relation [53]. It took several 
hours to compute Bfjj^ {I < 400), which is, however, a manageable one-time cost. 
We modified the CosmoMC package so that we can include the bispectrum likelihood (i.e. 
Eq. |4|) and parameters ({/^°£ al , U , 6 src }) in the Markov Chain Monte Carlo (MCMC) 
analysis [45, 54]. Note that the power spectrum likelihood is provided by the WMAP 
team and is already integrated into the CosmoMC package. In order to avoid unnecessary 
computations, we have added /nl an d 6 S rc a s amplitude parameters so that a radiation 



transfer function (i.e. gi(k) in Eq. |3.1| ) or the normalized bispectrum (i.e. Eq. 3J) and 



3.6| ) is not re-computed for a change of the parameters. For a cosmological model, we 
have considered ACDM + SZ effect + weak-lensing, where the cosmological parameters 
are A £ {f^, f2 c , r, n s , A s , A sz , Hq, f^ff 1 , /nl"'' ^src}- For data constraints, we have used the 
WMAP 7 year power spectrum and the bispectrum data described above. By running the 
modified CosmoMC package [45, 54] on an MPI cluster with 6 chains, we have explored the 
likelihood in the multi-dimensional parameter space. For a convergence criterion, we have 
adopted Gelman and Rubin's "variance of chain means" and set the R-l statistic to 0.03 for 
a stopping criterion [55,56]. Running a single CPU for each MPI process (total 6 CPUs), we 
reached the target convergence in less than a day. Once we reached the target convergence, 
we analyzed the MCMC samples with the CosmoMC package and obtained the parameter 
likelihood of /nl- In Fig. ||, we show the marginalized and mean likelihoods of /nl and b. 



src- 



The marginalized distribution shows the probability in the reduced dimension of parameter 
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Figure 4: Marginalized likelihoods in the plane of /nl versus cosmological parameters: the contour 
lines denote 1 and 2 cr intervals. 

space and the mean likelihood is associated with the best-fit value [54]. For a Gaussian 
distribution, marginalized and mean likelihoods are identical. A shown in Fig. there 
exist slight discrepancies between the marginalized and mean likelihoods, which indicate 
the deviation of the parameter likelihood from a Gaussian shape. In Fig. ||, we have 
plotted the marginalized likelihoods in the plane of /nl versus other parameters. As shown 
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Figure 5: MCMC samples: the x and y coordinates correspond to the values of /^l^ and /^l"' > 
color-coded according to the value of 6 src , marginalized over cosmological parameters. 



Table 1: The best-fit values and 68% confindence intervals for /nl and the cosmological parameters: 
(ACDM + sz + lens) 

parameter estimate 



flocal 

/NL oy -14 

U "76 < fZ U < * 

[KT 5 )uK 3 sr 2 ] -2-4±o9 

n b h* 0.0225t°oZl 

£l C h? O.H3^QQQg 

u - uoo -0.002 

n s 0.967i8;gu 

10 4 1 Q ioq+0.037 



logtlO 10 ^] 3.183_ 06 

A 1 4 o+0.58 

^*sz 1.42 

F [km/s/Mpc] 70.35±?;J 9 



in Fig. ||, there exist non-negligible deviations of the parameter likelihoods from a Gaussian 
shape. To be specific, the shape of parameter likelihood (e.g. contours of the same color) 
are poorly fitted by a Gaussian function. In Fig. ||, we show the MCMC samples, where 
the x and y coordinates correspond to the values of f^lf^ an< ^ /ni!" 1 ' anc ^ ^he va l ue °f ^src 
is color-coded. Note that MCMC samples shown in Fig. || are obtained after marginalized 
over the cosmological parameters. In Table [l], we show the best-fit values and la intervals 



- 11 - 



for the cosmological parameters as well as for /nl- As shown in Fig. [3| and Table [T], our 
result is in good agreement with the current results on /nl and the WMAP concordance 
model [22,28]. However, the confidence interval on f^f^ is different, which we attribute to 
a difference in the confidence interval estimation. While we have estimated the confidence 
interval by making a full exploration of the parameter likelihood, most of the existing works 
have relied on a Fisher matrix, which yields accurate results only for a parameter likelihood 
of a Gaussian shape. 

6. Discussion 

We have made an integrated MCMC analysis of the primordial non-Gaussianity (/nl)> us- 
ing the WMAP bispectrum and power spectrum. In our analysis, we have simultaneously 
constrained /nl and the cosmological parameters so that the uncertainties of the cosmo- 
logical parameters can properly propagate into the /nl estimation. From the parameter 
likelihoods deduced from the MCMC samples, we find that the parameter likelihood slightly 
deviates from a Gaussian shape, which makes the Fisher matrix estimation less accurate. 
Therefore, we have estimated the confidence intervals by exploring the parameter likelihood 
without using a Fisher matrix. We find that our best-fit values agree well with the existing 
results. However, we find that the confidence interval is slight different, which we attribute 
to the difference in the confidence estimation. Mostly due to the increase in the number of 
parameters, our approach requires a much higher computational load than the estimator 
of fixed cosmological parameters. However, with further optimization, our approach will 
be also feasible for the high angular resolution Planck data. 
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